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Abstract 

We present numerical studies of atomic transport in 3D and ID 
models for a mode-locked, pulsed atom laser as realized by Anderson 
and Kasevich [Science 281 (1998) 1686] using an elongated Bose con- 
densate of 87 Rb atoms poured into a vertical optical lattice. From our 
3D results we ascertain in a quantitative manner the role of mean- 
field interactions in determining the shape and the size of the pulses 
in the case of Gaussian transverse confinement. By comparison with 
ID simulations we single out a best-performing ID reduction of the 
mean- field interactions, which yields quantitatively useful predictions 
for all main features of the matter output. 
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1 Introduction 



The dynamics of Bose-Einstein condensates of alkali vapours 0-@ in very 
elongated traps is a matter of wide experimental and theoretical interest. 
An example is the mode-locked, pulsed atom laser which has been realized 
by Anderson and Kasevich Q by pouring a condensate of 87 Rb atoms in 
a vertical optical lattice. Drops of coherent matter leave the condensate 
under the effect of gravity. In this and other similar situations a theorist 
would like to reduce the problem to one-dimensional (ID) transport along the 
axial direction, keeping account of the interactions through a renormalization 
of the scattering length embodying the transverse confinement of the 3D 
condensate. 

A numerical study of atomic transport in a model relevant to the above- 
mentioned system has been based on the solution of a ID time-dependent 
Gross-Pitaevskii equation (GPE) [f|. This was obtained by freezing out the 
transverse motions of the condensate and by renormalizing the mean-field 
interactions according to a proposal made by Jackson et al. ||. The main 
result of this study was to show that, independently of the strength of the 
interactions and in complete agreement with the experimental data of An- 
derson and Kasevich M, the separation between successive matter drops is 
determined by the period of Bloch oscillations of the condensate in the ID pe- 
riodic potential of the optical lattice. It was further seen from the numerical 
results that the shape of the drops is closely related to that of the parent con- 
densate, thus supporting a mechanism of coherent emission and suggesting a 
practical way to tailor matter-wave laser pulses. Nevertheless, the question 
remains of how far a specific ID schematization may quantitatively capture 
other features of the phenomenon. 

In this Letter we address this question by carrying out a 3D numerical 
study of gravity-driven transport in a cylindrically symmetric condensate 
inside an optical lattice and by testing against these data the results of ID 
reductions of the problem. In addition to the proposal of Jackson et al. 0, we 
test a reduction previously proposed for condensates in anisotropic harmonic 
traps §, in which we fix the effective ID scattering length by adjusting the 
chemical potential of the ID model to that of the 3D condensate. As a further 
proof that the phenomenon of drop emission reflects the Bloch oscillations of 
the condensate, the spacing between the drops remains the same in all these 
models. The focus of our study then is on the shape and size of the matter 



2 



pulses. 

In our numerical simulations we use a fast, explicit time- marching scheme 
for the solution of the GPE in cylindrical geometry, which has been developed 
by Cerimele et al. 0. This method is briefly described in Section 2. Our 
extensive results of 3D simulation are presented in Section 3 and summarized 
in a simple diagram reporting the fractional number of particles in the first 
drop as a function of a single, suitably defined combination of system pa- 
rameters. The corresponding results from the ID reductions are reported in 
Section 4 and critically compared with those of the 3D simulation in Section 
5. This final section also present our conclusions. 



2 The numerical method 



The time-dependent GPE for the condensate wave function ^(r, t) is [|i0| , |TT 



where M is the atomic mass, Uj = A^h 2 aN/M is the interaction strength 
and U ex t is the external potential, a being the scattering length and N the 
number of particles in the condensate. In what follows U ext (r) is due to the 
optical lattice and to gravity, namely 

U ext (r, z) = U?[l - ex P (-r 2 /rf b ) cos 2 (2^/ A)] - Mgz , (2) 

Here, Uf is the well depth, rib is the transverse size, the wavelength A yields 
the lattice period d — A/2 and g is the acceleration of gravity. Finally, the 
normalization condition is / |^(r, t)| 2 <ir = 1. 

We numerically solve Eq. ([[]) by using an explicit time marching technique 
|J, in contrast with alternating-direction- implicit solver methods previously 
used in the context of Bose-Einstein condensation |T2], [Tj|. The present al- 



gorithm extends a fast time- staggered scheme proposed by Visscher |R| to 
solve the Schrodinger equation in an external potential, with the aims of 
preserving norm conservation in the presence of non-linear mean-field inter- 
actions and of handling cylindrical symmetry ||| . In brief, in solving Eq. 
([]]) we synchronously advance the real and imaginary parts of the scaled wave 
function in units of two time steps, using their intermediate centred value. 
The space derivatives are approximated by using centred differentiation. 
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The following crucial points deserve special comment. It is proven a priori 
||] and numerically verified that the algorithm preserves the norm of the wave 
function at each time step, provided that the boundary conditions are such 
as to annihilate surface terms. It can also be shown || that the numerical 
stability of the algorithm is preserved as long as the simulation time-step does 
not exceed a critical value Ar c , which is limited either by the grid spacing 
or by the magnitude of the product aN entering Uj. In our simulations the 
actual time-step is consistently kept well below the marginal stability thresh- 
old. Finally, the cylindrical symmetry is handled by an accurate treatment 
of the space derivatives near the symmetry axis [12| . 



Results for a 3D condensate in cylindrical 
geometry 



Equation (|I|) is made dimensionless by adopting the scale units Si = ^h/2Mui, 
S t = 1/uj and Se = fru for length, time and energy, ui being the radial fre- 
quency of the original magnetic trap. We also rescale the wave function by 
the dimensionless radial coordinate p. We use a grid resolution of 21 x 16 on 
each single well, the time step being At = 2 • 10~ 6 . 

The initial value \l/(r, z;t = 0) of the wave function is chosen so as to 
approximate the condensate realized in Q. Namely, 

exp[— Muj z [z — ld) 2 /2h] . (3) 

i 

Here, A is a normalization factor and I labels the occupied sites, their total 
number being n w . In constructing Eq. (Q) we have assumed (i) constant 
phase of the condensate in space [|15|; (ii) a Gaussian transverse profile with 



a frequency uj r = \J2Uf /(Mrf b ) given by the harmonic approximation to the 
transverse shape of the optical potential; and (iii) an overall axial profile 
reflecting that of the condensate inside the magnetic trap before loading the 
optical lattice and taken as a Gaussian having frequency Cj = 47r 3//5 /£ 2 C(j, 
with ( = (327rNa/ Si) 1 / 5 fT6f . Finally, the lowest state at each lattice site is 



occupied by a portion of condensate, having Gaussian shape with frequency 



ijj z = 2\JUf En/fi, where Er = h 2 /2M\ 2 is the recoil energy. 

Before entering the presentation of the simulation results, we list below 
the system parameters relevant to the experiment on 87 Rb H. These are 
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a = llOdo with cio the Bohr radius, N = 10 4 , A = 850 nm, = 80 /im, 
Uf = 1.4 Er and n w = 31. The parameters in this reference list provide our 
reference run. 

Fig. 1 shows four pictures of drop emission for different values of the 
coupling strength, by plotting the contour density profiles, taken with Uf = 
1.4 Eft after 5.3 ms. The first panel displays the behaviour of the non- 
interacting gas, namely the case a = and n w = 31. The other panels report 
the behaviour of the interacting gas with a = 110 a^; from left to right, the 
cases N = 10 4 with n w = 31, N = 10 5 with n w = 49 and N = 2 ■ 10 5 with 
n w = 57. All the other parameters are as listed above. 

A number of the results obtained in the earlier ID simulation H are 
recovered in the present 3D runs. Each drop in figure 1 extends over a number 
of wells equal to that occupied by the parent condensate. In all cases the 
drops are equally spaced by seventy wells from centre to centre. This spacing 
corresponds to 1.1 ms of simulation time, in agreement with experiment 
[fD and with the value of 1.09 ms for the period Tb = 2h/Mg\ of Bloch 
oscillations of the condensate in the periodic optical potential. The time 
lag between successive drops is independent of the amplitude of the periodic 
potential, of the strength of the interactions and of the dimensionality of 
the simulation sample, as expected if Bloch oscillations provide the correct 
interpretation of the observations |0|. Finally, both the axial width and 



the fine structure of each drop reproduce those of the parent condensate, as 
expected for coherent emission from all lattice wells. 

Again in analogy with the case of ID simulation [El, the transport be- 
haviour of the 3D system as a function of its governing parameters can be 
summarized in a single diagram. We introduce a scaling parameter g s , having 
the dimensions of the acceleration of gravity, by setting 

Mg s d = U? - Ui (4) 

where ^ = 47ih 2 aN/MR 3 with R = (S t is a measure of the mean-field inter- 
action strength. In figure 2 we plot the fractional number N^ rop /N of atoms 
in the first drop as a function of the ratio g/g s , each symbol representing 
a run with different input parameters. Symbols of different shape represent 
runs at different values of U^/Er, as shown in the legend. For each sym- 
bol we report the values of N drop /N corresponding to increasing coupling 
strength Ui, starting from the non- interacting gas and then increasing N in 
the sequence N = 10 4 , 10 5 and 2 • 10 5 (from left to right). 



5 



It is seen from figure 2 that there is a critical value for the onset of drop 
formation, which is g/g s ~ 0.14 (the onset is marked by an arrow). To all 
effects there is no emission of drops for subcritical values of g/g s , since the 
condition of resonance between the bound state in the well and the continuum 
is not satisfied. Well defined drops are instead emitted at supercritical values 
of g/g s . In this regime Nd rop /N increases rather regularly with g/g s and 
shows little sensitivity to the strength of the mean-field interactions at fixed 
Ui . Ultimately, with decreasing Uf the potential wells become too shallow 
and regular drop emission turns into a discharge of the whole condensate. 

On the other hand, an increase in the mean-field interactions affects the 
axial width of the drops. In particular the width of the first drop, as measured 
by its second moment, is close to n w /2, with an appreciable scatter from case 
to case, while the centre-to-centre distance of neighbouring drops is about 
70 wells. As a result, overlap between drops starts at iV ~ 3 • 10 5 . 

In summary, our 3D simulations confirm that the transport behaviour 
of a Bose condensate in a vertical optical lattice is described by a diagram 
as shown in figure 2. This diagram could yield useful predictions even for 
atomic species different from 87 Rb. 



4 One-dimensional modelling 

The time-dependent GPE for a ID reduction of the present transport problem 
is§ 

tk ^m^ = (-^+Ue X t(z)+U i m Z} t)\ 2) jij(z } t) , (5) 

where u ext {z) = Ui?,m 2, {2^z/\) — Mgz and ui = 4irh 2 aN/ M, a being a 
renormalized coupling parameter with the dimensions of an inverse length. 

In the renormalization proposed by Jackson et al. |J (hereafter referred 
to as model I) one assumes that the coherence length of the condensate in 
the axial direction is much larger than its transverse radius. The 3D wave 
function is factorized as \l/(r, t) = g(r, a(z{t))) where a(z) is the axial density. 
Using a harmonic approximation for the radial part of the optical potential, 
one obtains in the specific problem an effective scattering length a 1 = 7 7 a 



with 7 7 = yt/f I Er/ '(ribX). A similar renormalization of Uf is negligible, 
since ^> A. 
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As an alternative (hereafter referred to as model II) we impose that the 
value of the chemical potential of the 3D system be preserved upon reduction 
to ID [[7], ||. As is shown in a detailed calculation given in Appendix A 
within the Thomas-Fermi approximation, we find with 7 // given 

(for Uf > fi) by 

n 1 Uf-a 1 T . , 

~< = ^ + -Uf55D? J M • (6) 

an explicit expression for the positive quantity I(fi) being given in the Ap- 
pendix. We remark for later discussion that 7 7/ in Eq. @ depends implicitly 
on the product aN which determines the interaction strength. 

Table 1 reports the values of Sj and , y II Sf for a number of values of 
Uf /En and of N, all other parameters being as in the reference list. The 
values of h/Er are also shown. The important point of Table 1 is that while 
7 7 is independent of N, j 11 decreases with increasing N. This means that in 
model I an increase in the number of particles leads to a much more rapid 
increase of the mean-field interaction parameter. For instance, in the case 
Uf/Eji = 1.4 an increase of N by a factor 30 implies that , j II aN increases 
by only a factor 6.8. 

In the next Section we discuss the consequences of these different be- 
haviours of the two ID models. We compare the results of the 3D simulation 
reported in Section 3 with those obtained by solving the ID GPE with a 1 or 
a 11 inserted in turn into the mean-field term ui. 



5 Discussion and concluding remarks 

Table 2 collects the results for the fractional number Nd rop /N of atoms in 
the first drop for different values of Uf /E R and N in the interacting gas, as 
calculated from (i) the 3D simulation in cylindrical geometry (third column), 
(ii) the ID simulation in model I (fourth column), and (iii) the ID simulation 
in model II (last column). 

It is immediately seen that model II works better than model I and quan- 
titatively reproduces the data of the 3D simulation. The interactions lift the 
bound state towards the continuum by an amount which may be measured 
by the mean interaction energy Ej per particle . This is proportional to the 
product of the effective scattering length times the particle density. In the 
ID simulation according to model I we have Ej oc a/A oc a/A 2 r/fc. In the 3D 
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simulation we have instead Ej oc a/ (Xrf b ), which is significantly smaller since 
rib > A. 

Thus a picture emerges in which axial transport is accompanied by a 
transverse breathing of the condensate, due to the vanishing of radial con- 
finement at z = (2n + l)A/4 (see Eq. @). Model I cannot account for this 
behaviour, since it has been derived by freezing the transverse motion. It 
also follows that an increase in coupling strength is partially accomodated 
by a transverse spreading of the condensate, so that the value of Nd rop /N 
is rather insensitive to the mean-field interactions in the present range of 
system parameters. 

In summary, a Bose condensate in an optical lattice, subject to a con- 
stant driving field in the linear transport regime, behaves as a coherent blob 
of matter which executes Bloch oscillations through band states and can un- 
dergo Zener tunnel at the Brillouin zone edge. The simple diagram shown in 
Figure 2 describes the drops of coherent matter which are emitted via tunnel 
into the continuum, in dependence of the governing parameters of the sys- 
tem. The effect of the mean-field interactions is very small within the range 
of system parameters in which regular pulses are emitted. A quantitative 
reduction of this behaviour to a ID model can be achieved by imposing a 
simple condition of constancy of the chemical potential. We expect that this 
method of dimensionality reduction will be useful in other similar problems 
and applications. 
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A Calculation of the renormalization factor 

We set \I/(r, t) = exp(ifit/h)^/(r) in Eq. (|l|) and ip(z,t) = exp^ifit/fyiplz) in 
Eq. (|^). We then have to solve the stationary GPE's 

A**W = + U ext (r) + UjM^ *(r) (A.l) 

for \l/(r) in the 3D and 

^(z) = +u ext (z) + Ul mz)\^ ij(z) . (A.2) 

for ip(z) in ID. To this end we use the Thomas- Fermi approximation, which 
amounts to neglecting the kinetic energy terms. We then impose the nor- 
malization condition on \I/(r) and on ip(z), thereby obtaining the expressions 
for the chemical potential fi. 

After introducing dimensionless quantities as indicated in Section 3 and 
assuming confinement (Uf > yu), this procedure yields the relations 

32n2a ^ HS " = (1 _ 2(3) arccos (2/3 - 1) + y/l - (2(3 - l) 2 (A.3) 

in the ID case and 

32n 2 aNSf 32ir 2 aN-f n S? n „ rarcco S (2/3-i) 



/■arccos(/p — i) 

-2(3 / wtan(w/2)dw (A.4) 

Jo 



in the 3D case, where (3 = (Uf — n)/Uf. These equations yield Eq. @ in 
the main text, where 

I(u) = / wt&n(w/2)dw > (A.5) 
io 

with 

f(ji) = arccos(2([/° - //)/[/° - 1) . (A.6) 

In our numerical calculations we first determine \x from Eq. ( [A.l ) and then 
7 from Eq. (||). 
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Figure captions 

Fig. 1: Contour plots of the condensate density after 5.3 ms for Uf = 
1.4 Er, as functions of the axial coordinate z/d and of the transverse distance 
in fiin. The first panel on the left refers to the non-interacting gas with a 
number of particles iV = 10 4 . The other panels refer to the interacting gas, 
for a = 110 a and various values of N (from left to right, iV = 10 4 , 10 5 and 
2 • 10 5 ). 

Fig. 2: Diagram for drop formation from 3D simulation in cylindrical ge- 
ometry. We plot the fractional number N drop /N of particles in the first drop 
against g/g s for g — 981 cm/s 2 . The meaning of the symbols is explained in 
the legend and in the text. 
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Table 1: Calculated values of ^Sf, 7 // Sf and (J>/E R for various values of 
U?/E R and N. 
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Table 2: Fractional number Nd r0 p/N of atoms in the first drop for various 
values of Uf/E R and N. 
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